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ABSTRACT 

Molecular Cloud Complexes (MCCs) are highly structured and "turbulent" . Observational evidence 
suggests that MCCs are dynamically dominated systems, rather than quasi-equilibrium entities. The 
observed structure is more likely a consequence of the formation process rather than something that is 
imprinted after the formation of the MCC. Converging flows provide a natural mechanism to generate 
MCC structure. We present a detailed numerical analysis of this scenario. Our study addresses the 
evolution of a MCC from its birth in colliding atomic hydrogen flows up until the point when H2 
may begin to form. A combination of dynamical and thermal instabilities breaks up coherent flows 
efficiently, seeding the small-scale non-linear density perturbations necessary for local gravitational 
collapse and thus allowing (close to) instantaneous star formation. Many observed properties of MCCs 
come as a natural consequence of this formation scenario. Since converging flows are omnipresent in 
the ISM, we discuss the general applicability of this mechanism, from local star formation regions to 
galaxy mergers. 

Subject headings: turbulence — methods: numerical — ISM:molecular clouds 
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1. MOTIVATION 



Molecular Cloud Complexes (MCCs) in the Galaxy 
show an abundance of inte rnal structure, both in 
density and velocity (e.g. IStutzki fe Gustenl Il990t 
lFalgarondll99(l iMizuno et alJ 11995(1 . The density con- 
trasts are non-linear - independ ent of the importance 
of gravity l|Williams et alJll995l) -. and the MCCs ex- 
hibit non-thermal line- widths ^Falgarone fc Philipsil 990: 
Willi ams et alJ 1200(1 . generally interpreted as super- 
sonic turbulence. Filaments seem to dominate the 
morphologies in (column) density and velocity (e . 
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Bally et all 11987 



Churchwell et al 



Mizu no et all 119951 



2004) . The spatial distribution of den- 



sities and velocities is consistent with a turbulent spec- 
trum, the details of which are a matter of debate, though 
( Elm egreen fc Scaloll200l . 

Extremely puzzling is the source of this wealth of struc- 
ture. Wh ile stellar feedba ck certainly is a powerful driver 
(see e.g. IMac Low! 12004(1 it only acts locally, and defi- 
nitely only after the first stars have formed within the 
cloud. Moreover, stellar driving might be difficult to rec- 
oncile with the observed, nearly self-similar spatial en- 
ergy distribution. 

External drivers suffer from the fact that the cold 
dense gas esse ntially acts like a wall to any incom- 
mg wave (e.g. IVasauezI 119901 IBalsaralll996t lElmegreenl 
1999), preventing an efficient energy transfer from 
the warm diffu se component to the col d dense phase 
(see, however. iMiesch fc Zweibel 119941 for ID mod- 
els). iHennebelle fc I nutsukal 1)20051) argue that openings 
(holes, channels) in a MCC serve as an inroad for Alfvcn- 
waves, however, the question of what forms these chan- 
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nels remains. 

With growing evidence that the lifetimes of molec- 
ular clouds - at least in the solar neighborhood - 
range around a few ( 2 to 3) Myrs llElmegreenl I2000t 
lHartmann et all l200lt lHartmann! 12003). a picture in 
which MCCs are envisaged as transient objects in large- 
scale atomic colliding flows rather than as well-defined 
entities in a quasi-equilibrium state (Giant Molecular 
Cloud, GMC) is emerging. The colliding flows would 
accumulate atomic gas which might eventually reach col- 
umn densities high enough for ^-formation - at which 
point the "lifetime" of the molecular cloud would start, 
i.e. the accumulation time could be much longer. Some 



and form stars instantaneously 
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Elmegreenl 11993 


2000: 


Ballesteros-Paredes et^dj |1999 


Hartmann ct al. 


2001; 


Pringle et alJ 120011: lHartmann! 


2003). The large-scale 



flows might be driven by supernova explosions, Galactic 
shear, or might occur in galax y interactions (for r ecent 
evidence of cloud collisions, see lLoonev et al.l l2006). 

We propose and extend the idea (see that the struc- 
tures observed in MCCs arise from the very processes 
that form MCCs, i.e. they arise without recourse to non- 
linear perturbations put in by hand, an approach which 
defers the problem to an even earlier stage. We inves- 
tigate numerically the formation of structure in MCC 
precursors via colliding flows. To avoid the somewhat 
unwieldy "precursor of Molecular Cloud Complexes" , we 
will abbreviate this to PoMCloCs. By this we mean 
cold clouds (usually identified as Cold Neutral Medium, 
CNM), which may be fully atomic or may already con- 
tain traces of H2. In any case, the PoMCloCs serve as 
the initial stage of Molecular Cloud Complexes, MCCs. 

Our numerical models emphasize the ease with which 
colliding flows generate structure via a combination of 
dynamical and thermal instabilities. Due to the ther- 
mal instability, structure grows predominantly on small 
scales, leading to early non-linear density perturbations 
as possible seeds of gravitational collapse and star for- 
mation. The resulting line-widths in the cold gas when 
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seen in projection agree with observed values, however, 
the cold gas does not seem to exhibit internal supersonic 
turbulence. Only a few percent of the energy input re- 
mains available for driving turbulent motions in the cold 
gas, the bulk of the energy is lost due to radiative cool- 
ing. Within the scope of the models, we find that the 
time scale for H2 formation is limited not by a minimum 
temperature, but by the highly unsteady environment. 

Th is work extends our previous study ijHeitsch et al.1 
2005), in dimensions (here, we present corresponding 3D 
models), in resolution and in scope of physical applica- 
tions. Of course, we rely heavily on earlier work 
The physics of the problem is described in and the 
numerical realization as well as numerical artifacts are 
discussed in 21 an d ^6 6.11 in order to provide a back- 
ground for the interpretation of the results (SJSJ. These 
are summarized in while suggests possible routes 
to follow in the future. 

2. PREVIOUS AND CURRENT WORK 

The concept behind this study is that molecular cloud 
formation should be seen as a non-equilibrium process. 
When the perception of the ISM as a dynamical medium 
gained acceptance, many aspects of molecula r cloud and 
star formation t h eory w ere revisited. iHunterl l)1979|) and 
iHunter fc Fleck! l)1982D re-examined the Jeans criteria 
and found that converging velocity fields could provoke 
the gravitational co llapse of sub- Jeans mass gas clouds. 
Tohli ne~et al.l l)1987j) followed this line of thought, further 
emphasizing that diffuse gas cooling upon compression 
by even mild velocity and pressure disturbances could 
easily be transformed to high densities. They applied 
this idea to externally perturbed spherical clouds. Moti- 
vated by these investigations into the role of compressive 
velocity fields in promo t ing fa st cooling and gravitational 
collapse, Hun ter et al.l l)1986T) used two-dimensional nu- 
merical simulations to study the cooling and fragmenta- 
tion of gas compressed at the interface of two identical, 
oppositely directed supersonic colliding gas streams. 

Since these early works which explored some conse- 
quences of a dynamical ISM, molecular cloud forma- 
tion and evolution has been studied within the frame- 
work of compressible turbulence. A strongly fluctu- 
ating medium with motions on many scales results in 
high density albeit transient regions which might be 
identified as molecular clouds. Large-scale (kiloparsec) 
models of molecular cloud for mation were presented in 
two-dimensional simulatio ns bvlVazauez-Semadeni et alJ 
l)1995|) and iPassot et all l)1995|) . throwing light on the 
role of turbulence in the cycle of cloud and star forma- 
tion l)Larsonlll9811) . and specifying conditions for cloud 
formation in the presence of magnetic fields. But the nu- 
merical resolution needed to yield a detailed explanation 
of the structure and dynamics of molecular clouds makes 
it necessary to tackle the problem with smaller-scale sim- 
ulations (tens of parsecs). 

For instance, the role played by thermal instability con- 
tinues to be actively explored. This instability alone has 
been shown to be an efficient mechanism for generating: 
substructure and even driving turbulence ijBurkert fc Linl 
120001 iKritsuk & Normanll2f)02albL l2004|) As indicated 
by earlier work (e.s. iPriest fc Hevvaertslll974l in an ap- 
plication to the solar atmosphere), thermal instability 
can be triggered by compression. Indeed one-dimensional 



plane-parallel simulations of trans onic converging flows 
i|Hennebelle fc Peraultll 19991120001 ) show how this mech- 
anism transforms inflowing diffuse atomic gas to dense 
cold gas that is long-lived. Shock waves (driven 
e.g. by supernova explosions. iMac Low et alJ I20D51 
Ide Avillez fc BreitschwerdHl2005|) are another means to 
create cold, dense gas from warm, diffuse gas. Propagat- 
ing into the warm ISM, shock waves have been shown 
to fra gment in the presence of ther mal instability iFieldl 
[19651) and linear perturbations ijKovama fc Inutsukal 
I2OO0I I2002L I2004LD leading to H 2 -formation within a few 
Myrs l)Bergin et al.l 120041) . However the passage of a 
single shock might not leave in its wake enough dense, 
cold, cloudlets to constitute a MCC. Encounters between 
streams of transsonically transported gas, on the other 
hand might be a way to collect, compress, and cool large 
quantities of gas while at the same time possibly endow- 
ing the cold, dense gas with the dynamical and structural 
characteristics of molecular cloud complexes. This sce- 
nario has recently come into the spotlight. 

In high-resolution two-di mensional simulations of su- 
personic, colliding gas flows. lAudit fc Hennebelld ( 2005) 
confirm that even in a very dynamical environment a 
bistable medium develops as compressions initiate ther- 
mal instability. However their study concentrates on 
how turbulence generates and influences thermally un- 
stable gas, not on the reverse problem of how ther- 
mal instability could feed the turbulence characteristic 
of molecular clouds. In three-dimensional simulat i ons o f 
transonic colliding flows. IVazauez-Semadeni et alJ ( 2006) 
sh ow that a thin cold sheet, remi niscent of thos e observed 
by iHeiles fc Trolandl $200% and iHeilesI l|2004|) forms at 
the junction of the two flows. It develops turbulence 
which they a ttribute to the Non-linear Thin Shell Insta- 
bility (NTSI. lVishniacil9 941. Furthermore they find that 
even in simulations without gravity, the highest density 
gas is overpressured with respect to the mean warm neu- 
tral medium pressure, suggesting that the ram pressure 
of the colliding flows is responsible. 

To be fair, we should mention that favorable conditions 
for molecular cloud formation not only exist in converg- 
ing flows, but also in " focal p lanes" of (sustained) MHD- 
waves, as lElmegreenl (^999) shows. However, he notes 
that while highly structured MCCs form in such a sys- 
tem, the cold gas does not acquire significant turbulence. 

The notion of generating the turbulent substructure of 
molecular clouds by their form ation process has been dis- 
cussed in various conte xts (e.g. iBallesteros-Paredes et all 
1 19991 lHartmann et all 120011 lAudit fc Hennebelld 120051 
lyazouez-Semadeni et all I2006D. lAudit fc Hennebelld 
(2005) emphasized the evolution of the thermal states 
in the colliding flows, and provided a semi-analytical 
model for the evolution of a gas parcel. One of the 
main differences to the present study is that they im- 
pose vel ocity perturbations on th ei r inco ming gas flow 
- as do IVazauez-Semadeni et all (|2006), although as 
"random noise" at a percent-level - , whereas we de- 
fer th e structure generation to the ac tual flow collision 
site. IVazauez-Semadeni et all (120061) compared in de- 
tail analytical solution s with one-dimensiona l and t hree- 
dimensional models. iKovama fc Inutsukal l)2002l) and 
their subsequent work discuss molecular cloud formation 
behind a shock-compressed layer, more resembling the 
situation of an expanding shell. In their case, the initial 
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perturbations reside in the density field. The present 
study extends the "proof of concept" of "structure for- 
mation from scratc h" in two dimensions presented in 
iHeitsch et al.l l)2005fl by focusing on the dynamical state 
and the turbulent properties of the PoMCloC during its 
formation, and on the conditions necessary for the onset 
of H2 formation. 

3. PHYSICS 

We restrict the models to hydrodynamics with radia- 
tive cooling, leaving out the effects of gravity and mag- 
netic fields. Gravity would lead to further fragmentation, 
and magnetic fields could have a stabilizing effect (see 
also discussion in JFJ. For this regime, then, we iden- 
tify three relevant instabilities, namely the NTSI, the 
Kelvin-Hclmholtz-Instability (KHI) and the Thermal In- 
stability (TI) . This enumeration of course does not mean 
that the instabilities work independently of each other. 
Rather, the resulting dynamical system is a consequence 
of a combination of all three instabilities, however, in 
degrees depending on the local flow properties. 

3.1. NTSI 

The NTSI l)Vishniad ITool Fig. [TJ arises in a shock- 
bounded slab, when ripples in a two-dimensional slab fo- 
cus incoming shocked material and produce density fluc- 
tuations. Gas streaming along the cc-direction will be de- 
flected at perturbation peaks and collect in the troughs. 
Thus, x-momentum is transported laterally (along the y- 
direction). For a weak equation of state, the slab can act 
as a wall to deflect the incoming gas streams, while for 
an adiabatic one, the bounding shocks will travel faster 
at the y-locations of the troughs because of the higher 
com pression, thus leveling ou t the i nitial perturbations 
( see iVazauez-Se madcni et al.l (|2(306) for a discussion of 
the speed of the bounding shock). The growth rate is 
~ Csfc(fcA) 1 / 2 , where c s is the sound speed, k is the wave 
number along the slab, and A is the amplitude of the 
spatial perturbation of the slab. 

Various aspects of the NTSI h ave been investiga t ed nu- 
merically. iHunter et aD <|1986j) and iHueckstaedtl (|2003f) 
focused on the combination of the NTSI and self-gravity. 
The former authors found a criterion for forming self- 
gravitating fragments eit her under isotherm al condi- 
tions or with line-cooling. IHueckstaedtl ((2003) discussed 
the effect of the cooling strength on the evolution of 
the NTSI, stating that strong cooling leads to small- 
scale fragmentation whi ch wipes out the NTSI-effccts. 
iBlondin fe Mark's! l)1996j) investigated Vishniac's analy- 
sis numerically, noting that the shear flows along the slab 
lead to Kelvin-Helmholtz modes and identifying them as 
the main cause for the internal structur e of the slab. Al- 
thoug h in their radiative shock study, IWalder fc Folinil 
( 1998) focused mainly on the interaction of stellar winds, 
they extended the argument and the implications to 
the m o re general interstellar medium l)Walder fc Folinil 
IKlein fc Wooia lfT998h studied the NTSI in the 
context of cloud collisions, applying the process to the 
formation of clumpy filaments in e.g. the Orion Molec- 
ular Cloud. They partially included magnetic fields, fol- 
lowing the magnetic pressure term and the induction 
equation, but leaving out the tension term of the Lorentz- 
force. 



3.2. KHI 

The flows deflected at the slab (as indicated in Fig.^l 
can cause shear instability m odes, which have been stud- 
ied at great length (e.g. IChandrasekharl 119611: iSenl 
1961 iGerwinl Il963 iRobertsl 11974 iFerrari fc Trussonil 
19831 ISandham fc Reynolds! 119911 IKeopens et alJ 119991: 
Malagoli et alJll996l to name a few) . In the incompress- 
ible case, for a step function profile in the velocity and 
constant densities across the shear layer, the growth rate 
is given by the velocity difference kAU, i.e. all wave 
modes get unstable. If aligned with the flow, magnetic 
fields can stabilize against the KHI via the tension com- 
ponent of the Lorentz force, and the growth rate is given 
by k^J AU 2 — c A , i.e the modes that will get unstable are 
those for which the flow velocity difference is larger than 
the Alfven speed ca- For compressible gas, the situa- 
tion changes substantially: the system will be stable for 
all those wave numbers who se effective M ach number is 
larger than a critical value ljGerwinl il968'l. Thus, while 
Kelvin-Helmholtz modes definitely will play a role in the 
suggested scenario, Figure ^ certainly can give only a 
very simplified picture. 

3.3. TI 

The TI (jFieldl 1 1 965^ rests on the balancing of heating 
and cooling processes in the ISM. The TI develops an 
isobaric condensation mode and an acoustic mode, which 
- under ISM-conditions - is mostly damped. The con- 
densation mode's growth rate is independent of the wave 
length, however, since it is an isobaric mode, sma ller per- 
turbations will grow first l|Burkert fc Linll2000(l . Heat 
conduction sets the smallest growth scale. If this scale is 
not resolved nu merically, perturbations will grow at the 
resolution scale l|Kovama fc Inuts uka 2004a|). The signa- 
tures of the TI are fragmentation and clumping. They 
persist as long as the sound crossing time across a density 
perturbation is smaller than the cooling time scale 

3 kT 

|Tc| "2^AT (1) 

The TI can drive turbulence in an otherwise quiescent 
medium, ev en continuously, if an episodi c heating source 
is available (|Kritsuk fc Normadl2~002albl) . 

4. NUMERICS 

All three instabilities grow fastest or at least first 
on the smallest scales. This poses a dire chal- 
lenge for the numerical method. We chose a method 
base d on the 2nd order Bhatnagar-Gross-Krook formal- 
ism ilPrenderpst fc Xu1ll993t ISlvz fc Prendergastlll999t 
IHeitsch et al.ll2004USlvz et alJl2005j) . allowing control of 
viscosity and heat conduction. The code evolves the 
Navier-Stokes equations in their conservative form to sec- 
ond order in time and space. The hydrodynamical quan- 
tities are updated in time unsplit form. As shown below, 
statistical properties of the models are resolved with re- 
spect to grid resolution, viscosity and heat conduction, 
although the flow patterns change in detail — as is ex- 
pected in a turbulent environment. 

4.1. Heating and Cooling 

The heating and cooling rates are restricted to op - 
tically thin atomic lines following Wolfir e"et all {l995), 
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so that we are able to study the precursors of MCCs 
up to the point when they could form H2. Dust ex- 
tinction becomes important above column densities of 
iV(HI) w 1.2 x 10 21 cm~ 2 , which are reached only in the 
densest regions modeled. Thus, we use the unattenuate d 
UV radiation field for grain heating l|Wolfire et al] fl995h 
expecting substantial uncertainties in cooling rates only 
for the densest regions. The ionization degree is derived 
from a balance between ionization by cosmic rays and 
recombination, assuming that Ly a photons are directly 
reabsorbed. Numerically, heating and cooling is imple- 
mented iteratively as a source term for the internal en- 
ergy e of the form 



d t e = riT(T) - n 2 A(T) 



(2) 



in units of energy per volume per second. Here, T is 
the heating contribution (mainly photo-electric heating 
from grains), nA the cooling contribution (mainly due 
to the CII HFS line at 158/mi). Since the cooling and 
heating prescription has to be added outside the flux 
computations, it lowers the time order of the scheme. 
To speed up the calculations, equation J5J is tabulated 
on a 2048 2 grid in density and temperature. For each cell 
and iteration, the actual energy change is then bilinearly 
interpolated from this grid. 

4.2. Initial and Boundary Conditions 

Two opposing, uniform, identical flows in the x-y com- 
putational plane initially collide head-on at a sinusoidal 
interface with wave number k y = 1 (and k z = 1 for 3D 
models) and amplitude A (Fig. The incoming flows 
are in thermal equilibrium. The system is thermally un- 
stable for densities 1 < n < 10cm -3 . The cooling curve 
covers a density range of 10 -2 < n < 10 3 cm -3 and a 
temperature range of 30 < T < 1.8 x 10 4 K. The box 
side length is 44pc. Thus, Coriolis forces from Galactic 
rotation are negligible, however, the simulation domain is 
large enough to accommodate a moderately-sized molec- 
ular cloud complex. For an interface with A = 0, a cold 
high-density slab devoid of inner structure forms. The 
onset of the dynamical instabilities thus can be controlled 
by varying the amplitude of the interface perturbation. 
This allows us to test turbulence generation under mini- 
mally favorable conditions. This setup might seem artifi- 
cial: (1) The incoming flows are not expected to be per- 
fectly uniform, however, we chose to defer the moment 
of structure generation to the last possible moment in 
the simulation, instead of imposing perturbations on the 
incoming flows (see i)4 4.3|) . (2) As we will see below, the 
model run times extend considerably beyond lOMyrs. At 
an inflow speed of 10km s _1 this would correspond to a 
total extent of the system of 200pc. This suggests that 
the initial densities of the flows are more likely to be a 
few cm~ 3 to form molecular clouds with flows of order 
50-100 pc length. Spiral density waves can also produce 
coherent flows of the length required, at least in principle. 

The boundary conditions in the transverse (i.e. y and 
z) directions are periodic, while in the x-direction, the 
boundary values are defined as uniform inflow at con- 
stant density no and inflow speed vo. Consequently, the 
boundaries cannot treat outgoing waves or flows. Thus, 
the models have to be halted once the bounding shocks 
reach the boundaries. The corresponding time scale s 
have been discussed by Vaz auez-Semadeni et all |2006). 




X 

Fig. 1. — Sketch of the NTSI mechanism and geometry of the 
initial conditions. Gas is streaming along the x-direction and col- 
liding at the perturbed interface. The resulting shear flows excite 
KH-modes. 



To facilitate the analysis of the dynamics and time his- 
tory of the cold gas, the code is equipped with Lagrangian 
tracer particles, that are advected after each flux update. 
Tracer particles are deployed at grid cell centers where 
T < 300K and in four generations, starting at approxi- 
mately 5 Myrs in intervals of 2.5 Myrs. This allows us to 
study the effect of turbulent motions on the fluid parcels 
as well as the time history of (representat ive) regions of 
cold g as. Extending the model sequence of lHeitsch et alJ 
(2005), we varied the parameters, dimensions and reso- 
lution of the models (see Table . 

Before we discuss the simulation results and their phys- 
ical implications, we assess the numerical reliability of 
the code for the problem at hand. We will restrict the 
numerical discussion to the 2D models, since similar ef- 
fects pertain to the 3D models. 

4.3. Asymmetries 

FigureHshows stills of models 2C10c, 2C20c and 2C30c 
(i.e. a sequence in Mach numbers 1, 2 and 3) taken 11.5 
Myrs after flow contact time. While model 2C10c pre- 
serves the symmetry of the initial conditions, at larger 
Mach numbers, the structures develop strong asymmet- 
rical features with respect to the mid-plane. 

Since this could be cause for concern, we measured the 
degree of asymmetry in the models with time (Fig. |3J- 
Asymmetry is defined in terms of the density differences 
as 



An 



.1/2 



(3) 



where the average is over all cells with An > (and not 
over the whole domain in order to exclude the inflow ini- 
tial conditions which are symmetric), and An is the ab- 
solute value of the density difference between cells which 
should be symmetric across the upper and lower half of 
the simulation midplane. The initial conditions are per- 
fectly symmetric, but slight differences at the machine 
accuracy level in the reconstruction of hydrodynamical 
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Fig. 2. — Logarithmic temperature maps of the two-dimensional models 2C10c, 2C20c and 2C30c, i.e. with inflow speeds at Mach 1, 2 
and 3. The temperature range is given in logT[K], The stills are taken approximately 11.5 Myrs after flow contact. 



Table 1. Model parameters 



run 


N 


Mo 


vo [km s 1 ] 


To [K] 


no [cm 3 ] 


2B10 


b,c 


1.0 


5.3 


2.5 x 10 3 


1.0 


2B20 


b,c 


2.0 


10.6 


2.5 x 10 3 


1.0 


2B30 


b,c 


3.0 


15.9 


2.5 x 10 3 


1.0 


2C10 


b,c 


1.0 


8.9 


5.3 x 10 3 


0.5 


2C15 


b 


1.5 


13.4 


5.3 x 10 3 


0.5 


2C20 


b,c,d 


2.0 


17.8 


5.3 x 10 3 


0.5 


s/12C20 


b 


2.0 


17.8 


5.3 x 10 3 


0.5 


2C25 


b 


2.5 


22.4 


5.3 x 10 3 


0.5 


2C30 


b,c,d 


3.0 


26.7 


5.3 x 10 3 


0.5 


s/12C30 


b 


3.0 


26.7 


5.3 x 10 3 


0.5 


2C35 


b 


3.5 


31.2 


5.3 x 10 3 


0.5 


3C10 


a 


1.0 


8.9 


5.3 x 10 3 


0.5 


3C20 


a 


2.0 


17.8 


5.3 x 10 3 


0.5 


3C30 


a 


3.0 


26.7 


5.3 x 10 3 


0.5 



Note. — The first column lists the runs. The first digit of the 
model name stands for the dimension of the simulation, the second 
denotes the atomic particle density in the warm medium (where B 
corresponds to no = 1 cm -3 and C corresponds to no = 0.5 cm -3 
as indicated in the column marked no), and the third and fourth 
digits indicate the Mach number of the inflow. The second column 
lists the linear resolution N of the simulation where the letters a, 
b, c and d correspond to 256, 512, 1024 and 2048 grid cells per 
simulation box side. Where there is more than one letter, it means 
that more than one simulation was run (and not that the box size 
is oblong). The subsequent columns list the inflow Mach number 
.Mo, inflow velocity vo, initial temperature To and de nsity no. The 
models s/12C20 and s/12C30 are discussed in H66.ll 



variables eventually lead to a difference in the cooling 
strength and thus cause perceptible asymmetries. The 
code preserves perfect symmetry for a purely adiabatic 
equation of state (i.e. without the additional heating and 
cooling terms). After onset, the asymmetries grow lin- 
early in time until they reach a saturation level between 
A « 50 and 100. This corresponds approximately to the 
temperature (and density) contrast between the warm 
and cold gas and thus to the maximum asymmetry reach- 
able for the system. Although the asymmetries increase 
with Mach number and - to a lesser extent - with resolu- 
tion, they only appear well after the system has evolved . 
IKovama fe Inutsukal ll200l . lAudit k Hennebellel p005l) 
and lVazauez-Semadeni et alJ l)2006l) chose an alternative 
route: they added perturbations (in density or velocity) 



to the incoming flow, thus breaking the symmetry of the 
initial conditions. While the physical reason for adding 
perturbations to the inflow is perfectly obvious, in the 
present study we want to emphasize the point that even 
with the least possible perturbation to the flow, (tur- 
bulent) substructures are generated with ease. In some 
sense, this is an attempt to carry the argument to ex- 
tremes. 

Conversely, one could argue that imposing perturba- 
tions on the incoming flow helps to hide the symmetry 
breaking due to truncation errors. To get a stronger han- 
dle on how our somewhat extreme initial conditions are 
affecting the results, we repeated model 2C20 twice, once 
with an interface perturbation mode of k y — 32 instead 
of ky = 1, and in the second repetition with the k y = 32- 
mode overlaid on the k y = 1 mode at (l/32) 5 / 3 % of the 
amplitude at k y = 1, motivated by a turbulent cascade. 
In both cases, the system develops small-scale structures 
QuXi ky - — 32, meaning that the instability grows at the 
smallest imposed (k y — 32) scales, as long as these are 
larger than the Field length (see next section). 

4.4. Heat Conduction 

IKovama fc Inutsukal l)2004a|) demonstrated with ID 
models that the choice of the heat conductivity k can 
strongly influence the dynamics of the system. Espe- 
cially they argue that a too low heat conductivity can 
damp turbulent motions in the cold phase and can lead 
to arti ficial fra gmentation at t he grid scale (see however 
IVazauez-Semadeni et al.l l2006|) . The critical parameter 
here is the Field length l)Fieldlll965j) - the length scale 
below which heat conduction can stabilize the TI - 

^ 1/2 

v rt 2 A / 

where k is the heat conductivity (see below), and A 
is the cooling function from equation J2J. With a 
heat c onduction of k = 2 . 5 x 10 3 (T /[K]) 1 ' 2 l|Parkerl 
1953), IKovama fc Inutsukal l|2004af) conclude that for 
ISM-conditions, a linear resolution of several thousand 
cells is needed (they get convergence at 16384 cells). 

The heat conducti on in the BGK s cheme can be con- 
trolled explicitly fsee lSlvz et al.ll2006l for an analysis) by 
varying the kinematic viscosity, since the Prandtl num- 
ber in the code is Pr = 1 by construction. The choice of 
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against time for runs with resolution N = 512 (An; = 0.34pc). 



the viscosity is controlled by two considerations, namely 
(a) it should be large enough to prevent numerical arti- 
facts on small scales, and (b) it should be small enough 
to leave enough dynamical range. Obviously, both re- 
quirements are difficult to meet simultaneously. 

To establish to what degree our models are resolved, 
we begin by measuring the Field lengths (eq. 0]) cell- 
wise, excluding all cells which belong to the inflow, since 
their thermal time scale r c — > oo. Figure 0] shows the 
volume and mass fractions of the Field-unresolved cells 
against time for models 2C20c and 2C30c. Note that 
when the volume and mass fractions agree, the bulk of 
the non-equilibrium gas is in the cold phase (see also 
Fig. ll6|) . In other words, Figure0]and the following ones 
refer (mostly) to the cold gas phase. Thus, between 10% 
and 20% of the volume, and between 10% and up to 30% 
of the mass of the cold gas are not resolved. This is 
acceptable as long as we devise a selection criterion for 
the subsequent analysis. 

Figure El lets us estimate the size of the smallest struc- 
tures we can resolve. Again, we distinguish between vol- 
ume and mass fractions, which however in this case are 
close to identical. The figure indicates that the Field 
length for the bulk of the mass and volume lies between 
1 and 32 pc, scales which are well resolved by all our 
models. For our lowest resolution 2D runs (N = 512), 
the resolution limit is Ax — 0.34pc. Thus, for the subse- 
quent discussions, we will only consider structures with 
sizes larger than 0.34pc for 2D runs, and 0.68pc for 3D 
runs. 

If the Field-length is not resolved, structures tend to 
grow on grid scales, i.e with increasing resolution, there 



should be more and more small-scale structures. This ef- 
fect is clearly visible in the left panel of Figure El which 
shows a histogram of the size of cold regions (T < 300K) . 
The size of a cold region is defined as the geometric mean 
of its minimum and maximum diameter. The inset num- 
bers give the average number of cold regions per line of 
sight, i.e. normalized to resolution. If the Field length 
were resolved, heat conduction would lead to a cutoff at 
small L co id- Selecting for regions with L co u > 0.34pc, 
the histograms (and the average number of cold regions) 
agree sufficiently to proceed with the above selection cri- 
terion. 

5. RESULTS 

5.1. Turbulence 

Molecular Clouds consisten tly show non-thermal line- 
widths of a few km s -1 (e.g. iFalgarone 
Willi ams et alJl20uTih that - together with temperatures 
of T « 10K - are generally interpreted as supersonic 
turbulence. The line-widths in our models are consistent 
with the observed values (Fig. . The broad line wings 
are non-Gaussian. This may be a sign of inter mittency 
(e.g. IFalgarone fc PhillipsllT99l iLis et al.lll99(| . 

The "observed" lincwidth is derived from the his- 
togram of the density-weighted line-of-sight velocity dis- 
persion in the cold gas at T < 300K (Fig. El filled sym- 
bols). This linewidth would correspond to lincwidths 
in the cold neutral medium as e.g. traced by HI. 
Since the internal line-widths of coherent cold regions 
(open symbols) range around the sound speed of the 
cold gas (0.7km s" 1 ), the internal velocity disper- 
sions do not reach Mach numbers M. > 1 (see als o 
iKovama fc Inutsukal l200l lAudit fc Hennebelld IHU). 
Hence, the "supersonic" line- widths are a consequence 
of cold regions moving with respect to each other within 
a warmer and more diffuse medium, but not a result of 
internal supersonic turbulence in the col d gas which even- 
tually would be hosting star formation ijHartmann et af] 
l200lUKwan fc Sanderslll986ft . Note from Figure El that 
this result is independent of resolution and geometry. Be- 
cause of the thermal instability, to make figure El and El 
cold coherent regions are identified via a temperature 
threshold at T < 300K. Figure ITU1 shows the tempera- 
ture distribution according to volume and mass for mod- 
els 2C20 and 2C30 (at TV = 1024). Note that the above 
temperature criterion identifies the bulk of the cold gas. 
The choice of T = 300K as a "defining" temperature for 
the PoMCloCs is motivated by Figure Although the 
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Histograms of the Field lengths (eq. U) for models 2C20c (left) and 2C30c (right) 12 Myrs after flow contact 
denotes the resolution limit for N = 1024. 
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Fig. 6. — Histograms of the size of cold regions (T < 300K) per line of sight, i.e. normalized to resolution for models 2C20b, c, and 
d. The inset numbers give the average number of cold regions per line of sight. Left: All regions (down to the size of 1 cell) have been 
considered. There is a clear trend to larger numbers with increasing resolution. Right: Same histogram but selected for regions with 
Lcold > 0.34pc. 
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Fig. 7. — Histogram of the "observed" line-of-sight velocity dis- 
persion in the cold gas (T < 300K) for a model sequence in Mach 
numbers, 12 Myrs after initial flow contact. Shown are models at 
N = 512. Line-of-sight velocity dispersions are determined along 
the inflow direction (x-axis, see Fig. 0, and then laterally aver- 
aged. Observed velocity dispersion values are reached easily. See 
text for the choice of T = 300K as "defining" temperature for the 
PoMCloCs. 



histogram does not include gas at the inflow tempera- 
ture, the volume is dominated by warm gas (T ~ 10 3 ' 9 
K), while the mass is dominated by cold gas (T ~ 40 
K). The variations of the line-of-sight velocity dispersion 
with time is within the error bars shown. 

Are the line-widths actually indicating turbulent mo- 
tions, or are we seeing the inflow motions of (already) 
cold gas (e.g. iVazauez-Semadeni et al.ll200 6)? If the gas 
motions are truly turbulent, the average Lyapunov ex- 
ponents 



<A) = 



1 



t - 1, 



log 



d(t) 
d(to) 



(5) 



should be positive, indicating a growing separation be- 
tween neighboring particles. In equation to is the 
start time of particle advection, and d are the corre- 
sponding separations of neighboring particles. The av- 
erage refers to the simulation domain. After an initial 
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Fig. 8. — "Observed" (filled symbols) and internal (open sym- 
bols) velocity dispersion in the cold gas (T < 300K) against inflow 
speed for the same sequence of models as in Figure |7| 12 Myrs 
after initial flow contact. Resolution is N = 512. The dashed line 
denotes the sound speed in the cold gas. The internal velocity 
dispersion reaches Mach 1 at most. 
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Fig. 9. — "Observed" (filled symbols) and internal (open sym- 
bols) velocity dispersion in the cold gas (T < 300K) for a selection 
of models against inflow speed, 12 Myrs after initial flow contact. 
Question marks in the model names are wild cards for the Mach 
number. The dashed line denotes the sound speed in the cold gas. 
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Fig. 10. — Histogram of temperature for models 2C20 (left) and 2C30 (right) at N = 1024, and 12 Myrs after flow contact. The cold gas 
dominates the mass (thick lines), while the warm gas dominates the volume (thin lines). 



turbulent fFig. Illfl . With higher Mach number, the on- 
set of internal turbulence in the slab is delayed, which is 
mirrored by the lower Lyapunov exponents. Higher reso- 
lution leads to larger (A) at the time of the onset of tur- 
bulence, which is defined as t((A) > 0). For N = 1024, 
the initial separation of tracer particles is smaller, so that 
a larger spatial range can be covered. All exponents are 
positive, independent of resolution, and converge for late 
times. 

The power spectrum of the specific kinetic energy v 2 
(Fig. E| for model 3C20) can be used as another di- 
agnostic of turbulence in the models. Shown are the 
three linear spectra, taken along the inflow direction (a x ) 
and the transversal directions (a ViZ ). The spectral in- 
dex a x = —1.96 is consistent with the Fourier transform 
of a step function, as to be expected since the strong 
decelerations along the z-direction effectively lead to a 
discontinuity in v 2 (Note that for the initial conditions, 
not v 2 but v is discontinuous.). The spectral indices in 
the transversal directions, a y and a z , are consistent with 
a Kolmogorov spectrum, indicating fully developed tur- 
bulence. The lower spectrum (denoted by triangles and 
a corresponding slope a c — —3.24) shows the spectrum 
averaged over spheres of constant wavenumber k. Error 
bars denote errors on the mean. Within the errors, the 
slope is still consistent with a Kolmogorov spectrum of 
— 11/3 in three dimensions. The last half decade is dom- 
inated by numerical diffusion, leaving us approximately 
1.5 decades for physical analysis. 

The spectral indices differ depending on whether the 
spectrum is taken along the inflow or transversally. An- 
other way to see this is to split the specific kinetic en- 
ergy into a compressible (and translational) part with 
V x v = and a solenoidal part with V • v = 0. For 
the full domain (models 2C20 and 3C20), the specific 
kinetic energy is clearly dominated by the compressible 
part, because of the inflows (Fig.^J left). Note however, 
that the solenoidal part (thin lines) is steadily increasing. 
Even though the energy input rate is constant, the total 
specific kinetic energy (thick lines) decays with time as 
a consequence of the radiative cooling (see below). The 
specific energy restricted to the cold gas (at T < 300K, 
Fig. 1131 right) is dominated by the solenoidal component, 
and all components increase with time, indicating that 
once the gas has cooled down to its minimum tempera- 
ture, it starts to store the kinetic energy from the inflows. 
The solenoidal components are slightly larger for the 3D 
models, which is not surprising since the extra degree of 
freedom allows the gas to evade compression more effi- 
ciently. 

One might speculate as to whether the total specific 



kinetic energy will decrease until it has reached a mini- 
mum at the moment when (most of) the gas has cooled 
down to the minimum temperature given by the cool- 
ing curve, or whether it does find some kind of equilib- 
rium between kinetic energy input and thermal energy 
loss. Since with increasing density, the cooling will get 
more and more efficient, the latter seems unlikely. How- 
ever, it is probably even more unlikely that the inflows 
are maintained for times long enough to establish equi- 
librium between kinetic energy input and thermal en- 
ergy loss. Moreover, once gravity dominates the cold, 
dense regions, rather than reach a state of equilibrium, 
the clouds might pass through the phases of initial com- 
pression, turbulence generation, cooling and finally grav- 
itational collapse. 

5.2. Thermal states 

Structure forms in colliding flows as a result of an inter- 
play between dynamical and thermal instabilities. The 
dynamical instabilities generate compressions and shear 
flows, while the thermal instabilities amplify density per- 
turbations to the non-linear regime where they become 
potential seeds for self-gravitating structures. The effect 
of the TI on the thermal state of the gas is shown in 
Figure H3 

The solid line denotes the P(n) relation in thermal 
equilibrium. It is a direct result of the heating and 
cooling processes included and displays the usual three 
regimes, namely an "atomic" regime for logn < —0.5 
with an effective 7 e « 5/3, an "isothermal" regime for 
logn > 1.5 with j e — ► 1, and an unstable regime with 

7e < 0. 

Each cell in the scatterplots Figures HTI and [TBI is color- 
coded with its thermal time scale (eq. PP). Positive 
time scales denote heating (towards red colors) and neg- 
ative time scales correspond to cooling (towards blue col- 
ors). At high densities, the thermal time scales are short 
(greenish colors), which is mirrored in the very small 
scatter in log P around the equilibrium curve. Moving to 
lower densities, the cooling time scale increases dramati- 
cally, until we reach the instability regime, which consists 
of a black triangular region in P{n). This consists mostly 
of gas which has passed the bounding shocks (slightly 
increased densities due to the (nearly) adiabatic shock), 
but has not yet reached the cold dense phase. Since there 
is no gravity in the simulation, the ram pressure of the 
inflow (see dashed lines in Fig. I14fl limits the maximum 
pressure gas can reach in the system, apart from slight 
overshoots due to waves. A substantial amount of gas re- 
sides at higher pressure than that of the ambient inflow 
(higher by a factor of approximately 6 for model 2C20, 
and by 13 for model 2C30). This gas has been "overpres- 
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Fig. 11. — Left: Average Lyapunov exponents (eq. [F]) for a sequence of models in Mach number against time. After the initial 
compression phase, the cold gas becomes turbulent. Only the first tracer particle generation is used. Models shown in this diagram were 
run at N = 512. Right: Resolution effects. Higher resolution leads to larger (A), although all (A) > 0. 
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Fig. 12. — Specific kinetic energy spectrum for model 3C20 at 
15 Myrs. The linear spectral indices a IiSjZ refer to the coordinate 
directions, and a c gives the index of the spherically averaged spec- 
trum. Along the inflow, the strong decelerations effectively lead 
to a step function in < v 2 >, mirrored in the steeper index. Due 
to numerical dissipation, the last half decade is not available for 
physical analysis. 
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Fig. 13. — Left: Specific kinetic energy for the whole simulation 
domain of models 2C20 and 3C20 against time, split into compress- 
ible, solenoidal and total part. Compressible motions dominate be- 
cause of the inflows. Right: Specific kinetic energy for the cold gas 
(T < 300K) of models 2C20 and 3C20 against time, split into com- 
pressible solenoidal and total part. Solenoidal motions dominate 
in the cold gas. 



sured" by the highly compressible inflow. Even the gas in 
the isothermal branch (i.e. gas which has cooled down to 
the minimum temperature) resides at pressures between 
the initial thermal equilibr ium pressure and the ram pres- 
sure of the inflow (see also lVazauez-Semadeni et alJ200(il 
- an effect which traditionally h as been assigned to self- 
gravity, e.g. iPringle et al~l l2001). The more or less linear 
black regions at \ogn w —0.5 are a result of the adia- 
batic reaction of the shock-compressed gas, i.e. this gas 
has had not yet time to r eact to thermal effects (see also 
lYoung fe Mac Lowll2006j) . The dashed red line denotes 



the Hugoniot curve 

[ n 2 - 7 ■ - ni ) P 2 = [ m - 7 ; ] n 2 ) Pi, (6) 
V 7 + 1/ V 7+1 J 

that is, the pairs of values (P, n) for the state on one side 
(index 1) of the shock front that are compatible with the 
Rankine-Hu goniot jump conditions acro ss (index 2) the 
shock front llCourant fe Frie drichsl[l9 48). This has also 
been noted bv lAudit fc Hennebelld (|2005[) . 

From Figure El we see that the thermally unstable 
gas is generally still fast-moving, i.e. it has passed the 
bounding shocks but has not yet met the slab of dense, 
cold material (of course the initial inflow just shows up 
as a single dot in the P(n) and v(n) plots). High-density 
gas is mostly moving at a few km s _1 , consistent with 
the line-of-sight velocity dispersion discussed in i)5 5.11 

A small fraction of cells (1.7 x 10~ 3 for model 2C20, 
and 1.9 x 10~ 4 for model 2C30) seems to exhibit veloci- 
ties larger than the inflow velocity. This is hard to under- 
stand because there is no physical mechanism to acceler- 
ate gas to velocities higher than inflow velocities. All of 
those cells coincide with the largest density/temperature 
jumps occurring in the simulation and have slightly in- 
creased temperatures as well. Thus, they are a conse- 
quence of a slight overshoot from the reconstruction step 
in the scheme. Since the fraction of faulty cells does not 
vary over time (once the cold gas has formed) , we simply 
neglect these cells for the analysis. 

For our later discussion of critical masses for molecular 
cloud formation, Figure [TBI summarizes the mass content 
in the three thermal regimes, defined by the warm stable 
range (T > 3000K), the cold stable range (T < 300K) 
and the instability range in between. Obviously, once 
the cold regions form, most of the gas is locked there. 
The amount grows linearly with time, indicating that al- 
though the PoMCloC is turbulent, it essentially acts like 
a slab for collecting cold material. The warm phase (cen- 
ter panel) contains a close to constant amount of mass 
with time - it is just a transitory phase in our models. 
At any given time only a small fraction (but still a few 
percent) of the total mass is passing through the unsta- 
ble regime. This mass fraction depends slightly on the 
Mach number of the inflow, because the shorter dynam- 
ical time scales compete with the thermal time scales, so 
that gas which has been thrown out of thermal equilib- 
rium will have less chance of attaining equilibrium con- 
ditions again. A similar effect can be seen in the tem- 
perature maps (Fig. [5J: For lower Mach numbers, the 
transition from warm (yello w) to cold (blue) gas is much 
mor e pronounced (see also iSanchez-Salcedo et al.l 120021 
and iVazauez-Semadeni et alJl2006j) . Conversely, a high 
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Fig. 14. — Scatterplot of pressure against density 12 Myrs after initial flow contact, color-coded with the thermal time scale r c (eq.0, for 
models 2C20 and 2C30. Negative r c corresponds to cooling, positive t c to heating. The solid line denotes the equilibrium P(n)-distribution, 
the dotted lines stand for the initial conditions, and the dashed (black) line corresponds to the ram pressure of the inflow, p r = riv^. The 
dashed red line denotes the Hugoniot curve (eq. Irjl). 
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Fig. 15. — Scatterplot of absolute velocity against density 12 Myrs after initial flow contact, color-coded with the thermal time scale t c 
(eq.0, for models 2C20 and 2C30. Negative r c corresponds to cooling, positive r c to heating. Dotted lines denote the initial conditions. 
The fast-moving gas is generally in the thermally unstable regime. 



fracti on of gas in the unstable regime <|Heiles fc Trolandl 
12003^) might indicate a highly dynamical environment. A 
similar effect exists for 3D versus 2D models: Since the 
compressible modes are stronger in the 2D case, gas tends 
to be forced to cool in compressions, while in the 3D case, 
it can evade the compression and cooling by shearing mo- 
tions, which even would lead to additional heating. As 
Figure ITH1 shows, neither resolution nor geometry affects 
the mass fractions in the stable temperature regimes sig- 
nificantly. Especially, the fraction of mass in the cold 
phase (Fig. EH left panel) is indistinguishable for models 
at resolution N = 512 and N = 1024. 



5.3. Driving Efficiency 

Comparing the kinetic energy of the inflow to that in 
the unstable and cold (see i !5 5.2|) gas phases allows us 
to estimate the efficiency of turbulent "driving" in our 
models, i.e. the efficiency with which the highly ordered 
kinetic energy of the inflow is converted to the turbulent 
motions within the PoMCloC. Figure IT7I shows the ratio 
of the kinetic energy change to the energy input rate of 



the inflow, the efficiency, 



_ d t J nv 2 dV 
n v$A 



(7) 



in the cold (filled symbols) and unstable (open symbols) 
gas. The integral extends over all cells within the cho- 
sen temperature regime, and A is the area of the inflow. 
Overall, the efficiency £ in the cold gas ranges between 
2% and 5%. This is consistent with the ratio of the in- 
flow speed to the line-of-sight velocity dispersion (see e.g. 
Fig.^J. In other words, most of the energy is lost due to 
atomic line cooling. Gas in the unstable phase reaches 
an efficiency which is smaller by a factor of « 5: an effect 
of the lower densities (the velocities are generally higher, 
see Fig II 5t . There seems to be a weak trend to smaller 
efficiency £ with larger energy inflow uoVqA: Larger in- 
flow velocities result in higher compressions and thus in 
faster cooling, by which a growing fraction of the input 
energy is lost. 

5.4. Conditions for H2 Formation 
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Fig. 16.— Mass content in the warm (T > 3000K), unstable (300 < T < 3000K) and cold (T < 300K) phase, for a resolution (left panel) 
and geometry (right panel) test. Note that in order to compare between 2D and 3D models, the mass content is given per length, i.e. the 
3D models have been integrated along the ^-direction. For convenience, the mass fraction in the unstable regime has been multiplied by 
10. 
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Fig. 17. — Ratio of kinetic energy change to energy input rate in 
the unstable (open symbols) and cold (filled symbols) gas phases 
(see expression^ for the models indicated. Question marks in the 
model names are wild cards for the Mach numbers. Only a few 
percent of the energy input is transferred to the cold gas. 



A crucial point in our analysis is whether the cold gas 
reaches conditions favorable for H2 formation, and thus 
eventually for molecular clouds (so far we have only been 
talking about the precursors of molecular clouds). Three 
criteria need to be met. 

(1) For H2 formation, the gas temperature has to 
drop below a th reshold temperature which we take to 
be T c = 300K l|Cazaux et alJl200l . This might ap- 
pear to be somewhat high, since the sticking probabil- 
ity decreases rapidly with high er temperatures, and is 
near unity only fo r T < 20K ijCazaux fc Tielensl 120041 
iBergin et alJ 120041 . Nevertheless, we chose the above 
temperature limit, because it turns out that once the 
temperature of a gas parcel has dropped below T = 
300K, it will quickly cool down to the minimum temper- 
ature allowed by the cooling prescription. In that sense, 
the time scales give a lower limit, i.e. denote the point 
at which the first H2 could potentially be formed. Since 
the cooling curve only extends down to 30K anyway, we 
cannot make more detailed statements about H2 forma- 
tion. From Figure [TBI we already saw that for reasonable 
inflow speeds, a few 100M Q pc" 1 at T < 300K accumu- 
late within approximately 8 Myrs. This material would 
be in principle available for H2 formation - at least as 
far as the temperature is concerned. 

(2) The (now cold) gas has to stay cold "long" enough 
to allow for H2 formation. Time estimates for H2 for- 
mation vary. In their analysis of H2 formation behind 




Fig. 18. — Cumulative histogram of time intervals over which 
tracer particles stay at temperatures T < 300K (top), and over 
which tracer particles "see" column densities N > 10 21 cm -2 with 
respect to the ambient UV radiation field (bottom). Model se- 
quence in Mach number. Model 2C35 ended 6 Myrs after initial 
tracer deployment. All models were run at N = 512. 



shock fronts, IBergin et al.l l|2004|) quote time scales be- 
tween 5 and 10 Myrs, depending on the inflow momen- 
tum. Once H2 exists, its further formation could well be 
a run-away process, since self-shielding is much more ef- 
ficient than dust shielding. This means that the critical 
time scale is given by the onset of H2 formation. A min- 
imum requirement therefore is that the cold gas is not 
re-heated during this time scale. To measure this, we 
followed the temperature history of the tracer particles 
and determined how long each particle stays cold. The 
top panel of Figure ^] is a cumulative histogram of the 
time intervals over which tracer particles have temper- 
atures T < 300K. Apart from a small fraction at short 
time intervals AT (these are particles at the rims of the 
cold regions), most of the particles stay cold for at least 
6 Myrs. In fact, for most models, the cold gas parcels 
stay cold for over 14 Myrs, i.e. once it has cooled down, 
by far most of the gas stays cold. Model 2C35 ended 6 
Myrs after the initial tracer deployment, so in order to 
compare to the other models, all had to be evaluated to a 
maximum time interval of At — 6 Myrs, since otherwise 
it would look as if particles are re-heated in model 2C35 
after 6 Myrs. To summarize the top panel of Figure UhI 
the gas stays cold long enough to allow for H2 formation. 
(3) Finally, the cold, dense gas could be re-exposed to 
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the ambient UV radiation field. H2 formation requires a 
critical column density of N (HI) > 10 21 cm~ 2 . This we 
can determine again with the tracer particles, resulting in 
the bottom panel of Figure ITB1 which combines the tem- 
perature criterion of the top panel and the column den- 
sity threshold. Dropping the temperature criterion does 
not affect the result. Thus the critical quantity is the 
shielding column density, not the temperature. In other 
words, once the gas enters the "cold" phase, its thermal 
time scale is short compared to the dynamical time scale, 
so that the gas stays isothermal (see Fig. I14|) . However, 
due to the continuous re-structuring of the cloud, gas 
is repeatedly re-exposed to the UV radiation field and 
the column density of the cold gas only stays above 10 21 
cm~ 2 for about 1 Myr. This is a direct consequence of 
the PoMCloC's highly dynamical nature. Note that Fig- 
ure 1181 gives a pessimistic view: once a small fraction 
of the particles has reached conditions beneficial for H2 
formation, self-shielding will set in. The analysis up to 
now only includes shielding by dust. The inflow veloc- 
ity of model 2C10 is too low to reach sufficiently large 
column densities for H2 formation in the elapsed time 
of the simulation. As Berg in et alJ i|2004l) found in one 
dimension, average column densities must achieve values 
of 10 21 cm -2 for H2 formation to occur because of the 
shielding requirement. 

6. A SUMMARY AND DISCUSSION 

Molecular clouds in our Galaxy are complex and 
highly-structured, with broad, non-thermal line-widths 
suggesting substantial turbulent motions. Thus, molec- 
ular clouds very likely are not static entities and might 
not necessarily be in an equilibrium state, but their prop- 
erties could well be determined by their formation pro- 
cess. We presented numerical models of the formation 
of precursors of molecular cloud complexes (PoMCloCs), 
in large-scale colliding HI- flows. We will now summarize 
our findings and discuss their astrophysical relevance. 

6.1. Effects of Boundary Conditions 

In 21 we discussed to what extent numerical artifacts 
might influence our conclusions. One last effect needs to 
be mentioned, namely the choice of the boundary condi- 
tions. While the boundary conditions in the horizontal 
direction are prescribed by the inflow, we are free to de- 
fine the boundaries in the transverse direction. As men- 
tioned in i)44.2l our standard choice is periodic bound- 
ary conditions, i.e material leaving the simulation do- 
main at the bottom reenters at the top and vice versa 
(same holds for the tracer particles). This might raise 
the question of whether the level of turbulence and the 
amount of cold gas in the simulations is a "closed box" 
effect in the sense that if the boundaries were open the 
compressed gas could leave the simulation domain before 
cooling down and contributing to the cold gas mass. 

In order to assess the effects of this specific bound- 
ary choice, we ran a set of models with open bound- 
ary conditions in the direction transverse to the inflow, 
i.e. gas is free to leave the simulation domain. Mod- 
els s2C20 and s2C30 (Fig. top) have the same box 
size as models 2C20 and 2C30 (Fig. |2J, while in models 
12C20 and 12C30, the simulation domain was enlarged 
in the transversal direction (Fig. ^| bottom) in such a 
way that the inflow region has the same size as in the 




Fig. 19.— Top: Stills of models s2C20 and s2C30 with open 
boundary conditions in the transversal direction. The resolution 
is N = 512. Bottom: Stills of models 12C20 and 12C30 with open 
boundary conditions in the transversal direction and an "inactive" 
region above and below the inflow. The resolution is N = 512 X 
1024. 



square models (e.g.2C20, 2C30), but above and below 
this inflow region we place "inactive" regions into which 
the compressed gas can stream. In other words, not only 
are the boundaries in the transverse direction open, but 
there is no inflow in the x-direction in the "inactive" re- 
gions (in fact, they are open boundaries, i.e. material is 
free to leave the box). 

The morphologies of models with open and periodic 
boundaries are slightly different, but not disquietingly 
so. For the non-periodic models, note that while there is 
a clear outflow observable in 12C20 and 12C30, its signa- 
ture is weak in s2C20 and s2C30. The open boundaries 
are implemented via a constant extrapolation of the last 
active cells, so that the pressure gradient at the boundary 
in the transverse direction is constant. However because 
for models 12C20 and 12C30 we implement an "inactive" 
region above and below the inflow region, the gas in the 
dense compressed region, i.e. the inflow region, sees a 
pressure gradient along the transverse direction at the 
boundaries between the "active" and "inactive" regions. 

For models 12C20 and 12C30, dense gas is definitely 
squeezed between the colliding flows out into the "in- 
active" regions. How does this affect the mass budget? 
Figure E3 shows the mass per length contained in the 
three temperature regimes as in Figure IT6l but for mod- 
els 2C20, s2C20 and 12C20. 

Clearly, the total mass (i.e. the sum over the three 
temperature regimes) increases linearly, dominated by 
the growth of the cold mass fraction. Only model 12C20 
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Fig. 20. — Mass content in the warm (T > 3000K), unstable 
(300 < T < 3000K) and cold (T < 300K) phase, for models 2C20, 
s2C20 and 12C20. For convenience, the mass fraction in the un- 
stable regime has been multiplied by 10. The total mass in the 
simulation domain increases (more or less) linearly for all models, 
despite the open boundary conditions of models s2C20 and 12C20. 



shows some deviation from a linear evolution at later 
times. The cooling time scale is much shorter than the 
flow time scale and the sound crossing time, so that the 
gas is compressed and cooled down before it can feel the 
boundaries. In other words, gas in the incoming flow 
cools efficiently, and structures within the slabs form in- 
dependent of the boundary conditions. We conclude that 
our discussions of the mass budget and the level of tur- 
bulence are nearly independent of the specific boundary 
condition choice. 

6.2. Turbulence 

Observed line-of-sight velocity dispersions of a few km 
s _1 in the cold gas are reproduced by the models (Fig. 01 . 
Yet, the internal velocity dispersions are generally sub- 
sonic, i.e. it seems as if the term "supersonic" is not 
necessarily an accurate description of the hydrodynami- 
cal state of the cold gas (Figs.ElandEj!. IHartmannTl|2002fl 
argues, that because of the ages and small spatial disper- 
sions of young stars in Taurus, their velocity dispersions 
relative to their natal gas are very likely subsonic. The 
turbulent line-widths amount only to a fraction of the 
inflow velocity. 

The subsonic internal velocity dispersions come as a 
surprise, and the question arises whether self-gravity 
could lead to higher internal dispersions. While this 
is definitely something to test, to drive internal mo- 
tions, gravity would have to act on local non-linear 
perturbations within the cold gas, for whose growth 
there simply might not be enough time available 
( Burkcr tfc HartmannlkOO^ . 

The above result would affect some aspects of 
the model of turbulence-c ontrolled star formation (see 
iMac Low fc Klessenl 12004)) . There, supersonic turbu- 
lence leads to non-linear (isothermal) shock compressions 
which in turn trigger (at least partly) fragmentation and 
subsequent gravitational collapse. By this fragmentation 
mechanism star formation is "localized" in the sense that 
before the whole cloud can collapse, the non-linear den- 
sity perturbations caus ed by supersonic turbulence have 
formed stars (see also iBurkert fe Hartmannll2004|) . In 
the absence of supersonic turbulence in the cold regions, 
the non-linear density perturbations must arise during 
the formation process of the cloud. These perturbations 
could be envisaged as the cold filaments which are a con- 



sequence of the combination of dynamical and thermal 
instabilities. The subsonic internal motions are mirrored 
in the specific kinetic energy modes (Fig I13|) and spec- 
tra (Fig. I12|) . They clearly point to weakly compressible 
turbulence in the cold gas. 

Turbulence is generated with ease in the colliding flow 
scenario, and the resulting PoMCloC reproduces obser- 
vational bulk quantities. One possible consequence of 
turbulence as an initial-condition effect is the efficient 
mixing between warm and cold phases in the complex. 
There is an ongoing debate about the degree of mixing 
between atomic and molecular hydrogen in MCCs. While 
there is no doubt about the intermittent density (and 
velocity) structure in molecular clouds, it is less clear, 
what the diffuse, possibly warmer medium is made out 
of. Turbulent mixing would point to CNM, i.e. HI, while 
H2 self-shielding and H2 formation time scale arguments 
would favor a fully molecular cloud. The "openings" in 
the clouds (no matter whether they are e.g. channels, 
gaps or holes) might pla y a crucial role for the interna l 
dynamics of the clouds: iHennebelle fc Inutsukal l)2005fl 
argue that the channels of warm material in MCCs al- 
low energy (in their case in the form of Alfvcn waves) to 
enter the cloud and thus to drive the internal motions. 
However, this would only refer to motions of the cold gas 
regions with respect to each other, but not necessarily to 
internal motions in the cold gas. 

6.3. Thermal States 

Because of its very short thermal time scales, the dense 
gas closely follows the thermal equilibrium relation be- 
tween pressure and density (Fig. I14fl . This gas phase is 
approximately isothermal and corresponds to the dense, 
isolated regions which we identify as precursors of molec- 
ular gas. 

A perceptible amount of gas resides in the thcr- 
mally unstable phase, consistent with observations by 
iHeiles fc Trolandl |2003), also con sistent with the mod- 
eling results of other groups ( e.g. IDib fc BurkertJ 120051 
iVazauez-Semadeni et al.ll2006[) . The thermally unstable 
gas generally travels at high velocities (Fig. ll5|) . meaning 
that it corresponds to gas which has passed the bounding 
shock but has not yet "fallen onto" the cold, dense ob- 
jects (where it would cool down to the isothermal cloud 
temperature). The amount of unstable gas varies with 
the inflow Mach number, which is a consequence of com- 
peting dynamical and thermal time scales. 

6.4. Conditions for Hi-Formation 

Given a sufficiently high inflow speed, colliding flows 
are an efficient way to assemble enough cold gas to form 
molecular clouds. Despite the highly structured objects, 
the amount of cold mass increases approximately linearly 
with time as if collected in a cold slab (Fig. ll6|) . Once the 
gas has cooled below T « 300K in our models, it stays 
cold (Fig. I18|) : the cold material has very short ther- 
mal time scales, i.e. it is close to isothermal (in Fig. 
{d\ogP)/(d\ogn) -> 1 for n > 100 cm" 3 ). 

The main limiting factor for H2 formation in our mod- 
els is not the temperature of the gas, but the highly time- 
dependent column density. Because of the dynamical na- 
ture of the clouds, gas is continuously re-exposed to the 
ambient radiation field. Most of the gas manages to stay 
above shielding column densities of N > 10 21 cm -2 only 
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for 2Myrs at most (Fig. 1180. However, these estimates are 
pessimistic: (a) We only considered extinction by dust, 
not self-shielding by already formed H2. (b) Molecular 
clouds might not form "from scratch" , i.e. only from HI. 
Instead, the inflows mig ht already contain a substantial 
amount of H2 (see e.g. iPringle et al.l 120011) . This not 
only would shift the mass budget in favor of H2, but 
also would provide efficient self-shielding early on in the 
molecular cloud formation process. 

6.5. Cloud Formation Time Scales 

Even more than lOMyrs after flow contact time, in all 
our runs the total mass collected in the cold gas falls 
short of typical molecular cloud masses by a factor of 
approximately 10. The average column density essen- 
tially evolves one-dimensionally, i.e. N(t) — nv$t w 
1 x 10 20 cm _2 (i/Myr) for the fastest inflow speed in our 
models. Thus, at least 10 Myrs arc needed to reach 
an average column density high enough for self-shielding 
(since the dynamics lead to local focusing of the gas 
flows, isolated region s might reach this s tage earlier, but 
not markedly so, see lBergin et al.ll2004[) . Two remedies 
come to mind: (1) The cold mass growth rate depends 
nearly linearly on the inflow momentum (Fig. 116(1 . How- 
ever, the observed velocities and densities in the Galac- 
tic HI limits the inflow momentum to approximately 
no vq < 30km s _1 cm -3 , or - equivalently - no vq t s=s 10 21 
cm . Thus, molecular clouds might preferentially form 
in flows where the ambient atomic density is higher than 
the average density, i.e. no = 3 cm~ 3 , vq = 10 km s . 
Then, the mechanism described above might well domi- 
nate the formation of low-mass molecular clouds. More 
massive objects need either more time or a different pa- 
rameter regime as in e.g. galaxy mergers. (2) The models 
presented here do not account for gravity. Since a slab- 
like geometry is favored by construction in the colliding 
flow scenario, gravity would lead to collapse in the plane, 
causing the densities (and column densities) needed for 
shielding to be rapidly reached. In that sense, any time 
scales given for conditions favorable for H2 formation are 
upper estimates. Note that the lateral collapse can set in 
much later, and, if the cloud is large enough, and suffi- 
ciently dense substructures exist, might only occur after 
stars have formed locally. 

We would like to stress the point that since the pre- 
sented models are concerned with the formation of pre- 
cursors of molecular clouds, the somewhat largish as- 
sembly time scales do not consti tute a probl e m for 
short molecular cloud lifetimes. As IBergin "eTall ll200l 
discussed, the ratio of the precursor's assembly time 
scale over the actual molecular cloud lifetime can be as 
large as the relative amounts of atomic over molecular 
gas ( approximately a factor of 4 in t he solar neighbor- 
hood, lDamdll993l iSavage et al"lll977|) . For a molecular 
cloud life time of. say, 3Myrs in the solar neighborhood 
(|Hartmann et al.ll2001l and references therein) , the accu- 
mulation time for the P0MCI0C would be expected to 
range around > lOMyrs. Obviously, the molecular cloud 
lifetime starts only once molecules are being formed. 

Our models generally run longer than 10 Myrs. At an 
inflow velocity of 10 km s , this corresponds to a co- 
herent flow of approximately lOOpc, or strictly speaking, 
twice that length for a colliding flow. Although higher 
flow velocities and/or densities (note that at least with 



the densities we seem t o be on the lower s ide of the pa- 
rameter range, see e.g IBergin et alJ [20041) would allow 
shorter time scales and thus shorter coherence lengths, 
the question nevertheless arises whether such coherent 
flows are realistic. On smaller scales (and higher densi- 
ties), the interaction of winds o f nearby stars might pro - 
vide a possible source (see e.g. iChurchwell eTall 12004]) . 
while on larger scales, the most probable source would be 
interacting supernova shells. The somewhat longish as- 
sembly timescales suggest that the initial densities of the 
flows are more likely to be a few cm~ 3 to form molecular 
clouds with flows of order 50-100 pc length. Spiral den- 
sity waves can also produce coherent flows of the required 
length, at least in principle. 

7. FUTURE WORK 

The main purpose of this work was to demonstrate 
that precursors of molecular cloud complexes forming in 
colliding flows reproduce observational measures such as 
line-widths and column densities, without any recourse 
to initially imposed perturbations. The combination of 
dynamical and thermal instabilities efficiently generates 
the non-linear substructures necessary for a fast gravita- 
tional collapse and a (close to) instantaneous onset of star 
formation, thus meet ing the short tim e scales required 
by observations (e.g. iHartm ami 2003). Once the mate- 
rial is cold and dense, it seems difficult to drive highly 
supersonic moti ons (unless there is an intern al driver, 
mea ning stars') (iVascmezI 119901: lBalsaralll996L see how- 
ever |Miesch^^Z^ibeM9^Jor even to mix it with warm 
gas ((Heitsch et al. 2006). However, the cold and warm 
gas must be efficiently mixed in order to allow the cold 
material to move "freely" at velocities which are super- 
sonic with respect to the cold gas' temperature. Thus, 
the cold substructures must arise during the formation 
process of the cloud, from compression and cooling in 
the warm phase. We demonstrated that colliding flows 
provide a generic physical mechanism for creating non- 
linear structure in PoMCloCs, without any recourse to 
perturbed initial conditions. In their further evolution 
(i.e beyond the models presented here), these non-linear 
perturbations will be the seeds for gravitational collapse 
and star formation. Thus, our models show that - from 
a dynamical point of view - the concept of "short" cloud 
life times is feasible. Ultimately, one can begin to see how 
realistic the concept of "turb ulent support" of molecu- 
lar cl ouds actually is (see e.g. iBallesteros-Paredes et alJ 
1999). 

An MCC in the solar neighborhood has an aver age co l- 
umn density of N(BI) « 1.4 x 10 21 cm" 2 iMcKedll999|) . 
coinciding with the dust-shielding column density neces- 
sary for H2 formation. Once this value is reached, H2- 
formation is expected to begin. Within 2 or 3 more Myrs, 
star formation will have occurred, and the cloud will have 
dispersed again. Within this scenario, the average col- 
umn density of an ensemble of MCCs should be roughly 
constant, implying a mass-size relation of M oc L 2 . The 
latter would be expected for a lLarsorJ l|1981|) relation of 
the type a oc L 2 / 5 , c onsistent with observed scalings (e.g. 
Elm egreen fc Scaldl2004|) . 

Clearly, this study leaves us with many unanswered 
questions. Gravity might lead to global edge effects 
(Burk ert fc Hartmannl 12004) as well as to local "super- 
sonic" velocities and to additional instabilities (see e.g. 
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iHunter et aDll997t iHueckstaedt et alJl2005j) . Actual H 2 
formation would further lower the temperature in the 
cold gas, in which case eventually supersonic velocities 
might be reached. The inflows might already be partially 
molecul ar, in which case t he filling factor of H2 should 
be low (Pringl e~et aUl200lt) . And finally, magnetic fields 
could influ ence and maybe control the structure forma- 
tion (e.g. iVazquez-Semadeni et al.l Il995t iPassot et al.1 
ll995UElmegreenlll999HHartmann et al]l2001ft . ~ 
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